---
title: "Casas-figure-9"
format: html
editor: visual
---

```{r}
library(devtools)
install.packages("lme4", type = "source")
install_github("naoki-egami/dsl", ref = "update-Aug2024", dependencies = TRUE, force = TRUE)
library(dsl)
library(janitor)
library(tidyverse)
library(broom)
library(jtools)
library(tidyverse)
library(ggplot2)
library(estimatr)
library(broom.mixed)
```

## Making Figure 9 (Casas, $N = 1200$, $k=500$)

```{r}
gpt4o <- read_csv("~/Downloads/DSL images/gpt4o.csv")
itm <- read_csv("~/Downloads/DSL images/itm-or-tweets-w-img-fname.csv")
```

```{r}
itm <- itm %>% filter(!is.na(img_fname)) %>% filter(img_fname %in% unique(gpt4o$name))
length(unique(itm$img_fname))
nrow(drop_na(gpt4o))
gpt4o <- drop_na(gpt4o)
```

```{r}
casas_sampling <- function(itm, gpt4, N, seed){
  set.seed(seed)
  casas_merged_0 <- itm %>% 
  left_join(gpt4, by = c("img_fname" = "name"))
  
  casas_merged_0[,c("anger.x", "fear.x", "disgust.x", "sadness.x", "enthusiasm.x")] <- casas_merged_0[,c("anger.x", "fear.x", "disgust.x", "sadness.x", "enthusiasm.x")]/10
  casas_merged_0[,c("anger.y", "fear.y", "disgust.y", "sadness.y", "enthusiasm.y")] <- casas_merged_0[,c("anger.y", "fear.y", "disgust.y", "sadness.y", "enthusiasm.y")]/10


  # Sample with replacement from the rows of the data frame
  sample_casas <- casas_merged_0[sample(nrow(casas_merged_0), size = nrow(casas_merged_0), replace = TRUE), ]

  # Randomly select a subset of rows to assign NA (exclude N rows)
  sub_sample_indices <- sample(nrow(sample_casas), size = (nrow(sample_casas) - N))
  sub_sample_casas <- sample_casas
  sub_sample_casas$row <- seq_len(nrow(sub_sample_casas))

  # Assign NA to our labeled observations 
  sub_sample_casas$anger.x[sub_sample_indices] <- NA
  sub_sample_casas$fear.x[sub_sample_indices] <- NA
  sub_sample_casas$disgust.x[sub_sample_indices] <- NA
  sub_sample_casas$sadness.x[sub_sample_indices] <- NA
  sub_sample_casas$enthusiasm.x[sub_sample_indices] <- NA
  return(sub_sample_casas)
}
```

```{r}
# number of times we want to run this
k = 500
# number of observations we want to sample with replacement
N = 1400

# set seed for replicability
set.seed(23928384)

seeds <- sample(x=1:99999,size=k)

```

```{r}


# make data frames where we will save our models 

dsl_anger_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
dsl_disgust_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
dsl_enthusiasm_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
dsl_fear_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
dsl_protest_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
dsl_sadness_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))
dsl_symbol_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))


sub_anger_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
sub_disgust_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
sub_enthusiasm_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
sub_fear_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
sub_protest_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
sub_sadness_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))
sub_symbol_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))

gpt_anger_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
gpt_disgust_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
gpt_enthusiasm_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
gpt_fear_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
gpt_protest_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate", "std.error", "conf.low", "conf.high", "p.value"))
gpt_sadness_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))
gpt_symbol_1200 <- data.frame(matrix(nrow = 5,ncol = k), row.names = c("estimate","std.error", "conf.low", "conf.high", "p.value"))
```

```{r}
for (i in 1:k) {
  # Generate the data sample
  df <- casas_sampling(itm, gpt4o, N, seed = seeds[i])
  
  # Initialize variables for DSL results
  d <- NULL
  
  # Wrap DSL model fitting and saving results in a tryCatch block
  tryCatch({
    dsl_model <- dsl(
      model = "lm", 
      formula = log(retweet_n+1) ~ followers_count +
        friends_count +
        prev_tweets + time_control +
        protest +
        symbol +
        anger.x +
        fear.x +
        disgust.x +
        sadness.x +
        enthusiasm.x,
      predicted_var = c("anger.x", "fear.x", "disgust.x", "sadness.x", "enthusiasm.x"),
      prediction = c("anger.y", "fear.y", "disgust.y", "sadness.y", "enthusiasm.y"),
      data = df, tuning = TRUE
    )
    
    # Save DSL results
    d <- data.frame(summary(dsl_model))
    dsl_protest_1200[,i] <- c(d$Estimate[10], d$Std..Error[10], d$CI.Lower[10], d$CI.Upper[10], d$p.value[10])
    dsl_symbol_1200[,i] <- c(d$Estimate[11], d$Std..Error[11], d$CI.Lower[11], d$CI.Upper[11], d$p.value[11])
    dsl_anger_1200[,i] <- c(d$Estimate[12], d$Std..Error[12], d$CI.Lower[12], d$CI.Upper[12], d$p.value[12])
    dsl_fear_1200[,i] <- c(d$Estimate[13], d$Std..Error[13], d$CI.Lower[13], d$CI.Upper[13], d$p.value[13])
    dsl_disgust_1200[,i] <- c(d$Estimate[14], d$Std..Error[14], d$CI.Lower[14], d$CI.Upper[14], d$p.value[14])
    dsl_sadness_1200[,i] <- c(d$Estimate[15], d$Std..Error[15], d$CI.Lower[15], d$CI.Upper[15], d$p.value[15])
    dsl_enthusiasm_1200[,i] <- c(d$Estimate[16], d$Std..Error[16], d$CI.Lower[16], d$CI.Upper[16], d$p.value[16])
  }, error = function(e) {
    message(sprintf("Error in DSL model at iteration %d: %s", i, e$message))
  })
  
  # Reset DSL variables
  dsl_model <- NA
  d <- NA
  
  # Fit subsample human model
  subsample_human_model <- lm_robust(log(retweet_n+1) ~ followers_count +
                                 friends_count +
                                 prev_tweets +
                                 time_control +
                                 protest +
                                 symbol +
                                 anger.x +
                                 fear.x +
                                 disgust.x +
                                 sadness.x +
                                 enthusiasm.x, data = df)
  
  d <- subsample_human_model %>% tidy(conf.int = TRUE)
  d <- as.data.frame(d)
  
  # Save results from the subsample human model
  sub_protest_1200[,i] <- c(d$estimate[10], d$std.error[10], d$conf.low[10], d$conf.high[10], d$p.value[10])
  sub_symbol_1200[,i] <- c(d$estimate[11], d$std.error[11], d$conf.low[11], d$conf.high[11], d$p.value[11])
  sub_anger_1200[,i] <- c(d$estimate[12], d$std.error[12], d$conf.low[12], d$conf.high[12], d$p.value[12])
  sub_fear_1200[,i] <- c(d$estimate[13], d$std.error[13], d$conf.low[13], d$conf.high[13], d$p.value[13])
  sub_disgust_1200[,i] <- c(d$estimate[14], d$std.error[14], d$conf.low[14], d$conf.high[14], d$p.value[14])
  sub_sadness_1200[,i] <- c(d$estimate[15], d$std.error[15], d$conf.low[15], d$conf.high[15], d$p.value[15])
  sub_enthusiasm_1200[,i] <- c(d$estimate[16], d$std.error[16], d$conf.low[16], d$conf.high[16], d$p.value[16])
  
  # Fit GPT4 model
  gpt4_model <- lm_robust(log(retweet_n+1) ~ followers_count +
                      friends_count +
                      prev_tweets +
                      time_control +
                      protest +
                      symbol +
                      anger.y +
                      fear.y +
                      disgust.y +
                      sadness.y +
                      enthusiasm.y, data = df)
  
  d <- gpt4_model %>% tidy(conf.int = TRUE)
  d <- as.data.frame(d)
  
  # Save results from the GPT4 model
  gpt_protest_1200[,i] <- c(d$estimate[10], d$std.error[10], d$conf.low[10], d$conf.high[10], d$p.value[10])
  gpt_symbol_1200[,i] <- c(d$estimate[11], d$std.error[11], d$conf.low[11], d$conf.high[11], d$p.value[11])
  gpt_anger_1200[,i] <- c(d$estimate[12], d$std.error[12], d$conf.low[12], d$conf.high[12], d$p.value[12])
  gpt_fear_1200[,i] <- c(d$estimate[13], d$std.error[13], d$conf.low[13], d$conf.high[13], d$p.value[13])
  gpt_disgust_1200[,i] <- c(d$estimate[14], d$std.error[14], d$conf.low[14], d$conf.high[14], d$p.value[14])
  gpt_sadness_1200[,i] <- c(d$estimate[15], d$std.error[15], d$conf.low[15], d$conf.high[15], d$p.value[15])
  gpt_enthusiasm_1200[,i] <- c(d$estimate[16], d$std.error[16], d$conf.low[16], d$conf.high[16], d$p.value[16])
  
  d <- NA
  print(i)
}

```

## Saving results and making plots

```{r}

dsl_protest <- as.data.frame(t(dsl_protest_1200))
dsl_symbol <- as.data.frame(t(dsl_symbol_1200))
dsl_anger <- as.data.frame(t(dsl_anger_1200))
dsl_fear <- as.data.frame(t(dsl_fear_1200))
dsl_disgust <- as.data.frame(t(dsl_disgust_1200))
dsl_sadness <- as.data.frame(t(dsl_sadness_1200))
dsl_enthusiasm <- as.data.frame(t(dsl_enthusiasm_1200))

dsl_protest  <- na.omit(dsl_protest)
dsl_symbol  <- na.omit(dsl_symbol)
dsl_anger  <- na.omit(dsl_anger)
dsl_fear  <- na.omit(dsl_fear)
dsl_disgust  <- na.omit(dsl_disgust)
dsl_sadness  <- na.omit(dsl_sadness)
dsl_enthusiasm  <- na.omit(dsl_enthusiasm)



## get Benchmark estimates

org <- itm  
org[,c("anger", "fear", "disgust", "sadness", "enthusiasm")] <- org[,c("anger", "fear", "disgust", "sadness", "enthusiasm")]/10

human_model <- lm_robust(log(retweet_n+1) ~ followers_count +
                              friends_count +
                              prev_tweets +
                              time_control +
                              protest +
                              symbol +
                              anger +
                              fear +
                              disgust +
                              sadness +
                              enthusiasm, data = org)

h_tidy <- human_model  %>% tidy(conf.int = TRUE)

b <- as.data.frame(h_tidy) 


t_protest <- b$estimate[10]
t_symbol <- b$estimate[11]
t_anger <- b$estimate[12]
t_fear <- b$estimate[13]
t_disgust <- b$estimate[14]
t_sadness <- b$estimate[15]
t_enthusiasm <- b$estimate[16]


# calculate percentage covered

dsl_protest["covered"] <-  with(dsl_protest, t_protest >= conf.low & t_protest <= conf.high)
dsl_symbol["covered"] <- with(dsl_symbol , t_symbol >= conf.low & t_symbol <= conf.high)
dsl_anger["covered"] <- with(dsl_anger , t_anger >= conf.low & t_anger <= conf.high)
dsl_fear["covered"] <- with(dsl_fear , t_fear >= conf.low & t_fear <= conf.high)
dsl_disgust["covered"] <- with(dsl_disgust , t_disgust >= conf.low & t_disgust <= conf.high)
dsl_sadness["covered"] <-with(dsl_sadness , t_sadness >= conf.low & t_sadness <= conf.high)
dsl_enthusiasm["covered"] <-with(dsl_enthusiasm, t_enthusiasm >= conf.low & t_enthusiasm <= conf.high)


percent_dsl_protest  <- sum(dsl_protest["covered"], na.rm = TRUE)/nrow(dsl_protest["covered"])
percent_dsl_symbol  <- sum(dsl_symbol["covered"], na.rm = TRUE)/nrow(dsl_symbol["covered"])
percent_dsl_anger  <- sum(dsl_anger["covered"], na.rm = TRUE)/nrow(dsl_anger["covered"])
percent_dsl_fear  <- sum(dsl_fear["covered"], na.rm = TRUE)/nrow(dsl_fear["covered"])
percent_dsl_disgust  <- sum(dsl_disgust["covered"], na.rm = TRUE)/nrow(dsl_disgust["covered"])
percent_dsl_sadness  <- sum(dsl_sadness["covered"], na.rm = TRUE)/nrow(dsl_sadness["covered"])
percent_dsl_enthusiasm  <- sum(dsl_enthusiasm["covered"], na.rm = TRUE)/nrow(dsl_enthusiasm["covered"])




# coverage for GPT

gpt_protest <- as.data.frame(t(gpt_protest_1200))
gpt_symbol <- as.data.frame(t(gpt_symbol_1200))
gpt_anger <- as.data.frame(t(gpt_anger_1200))
gpt_fear <- as.data.frame(t(gpt_fear_1200))
gpt_disgust <- as.data.frame(t(gpt_disgust_1200))
gpt_sadness <- as.data.frame(t(gpt_sadness_1200))
gpt_enthusiasm <- as.data.frame(t(gpt_enthusiasm_1200))

gpt_protest  <- na.omit(gpt_protest)
gpt_symbol  <- na.omit(gpt_symbol)
gpt_anger  <- na.omit(gpt_anger)
gpt_fear  <- na.omit(gpt_fear)
gpt_disgust  <- na.omit(gpt_disgust)
gpt_sadness  <- na.omit(gpt_sadness)
gpt_enthusiasm  <- na.omit(gpt_enthusiasm)



gpt_protest["covered"] <-  with(gpt_protest, t_protest >= conf.low & t_protest <= conf.high)
gpt_symbol["covered"] <- with(gpt_symbol, t_symbol >= conf.low & t_symbol <= conf.high)
gpt_anger["covered"] <- with(gpt_anger, t_anger >= conf.low & t_anger <= conf.high)
gpt_fear["covered"] <- with(gpt_fear, t_fear >= conf.low & t_fear <= conf.high)
gpt_disgust["covered"] <- with(gpt_disgust, t_disgust >= conf.low & t_disgust <= conf.high)
gpt_sadness["covered"] <-with(gpt_sadness, t_sadness >= conf.low & t_sadness <= conf.high)
gpt_enthusiasm["covered"] <-with(gpt_enthusiasm, t_enthusiasm >= conf.low & t_enthusiasm <= conf.high)

percent_m_protest <- sum(gpt_protest["covered"], na.rm = TRUE)/nrow(gpt_protest["covered"])
percent_m_symbol <- sum(gpt_symbol["covered"], na.rm = TRUE)/nrow(gpt_symbol["covered"])
percent_m_anger <- sum(gpt_anger["covered"], na.rm = TRUE)/nrow(gpt_anger["covered"])
percent_m_fear <- sum(gpt_fear["covered"], na.rm = TRUE)/nrow(gpt_fear["covered"])
percent_m_disgust <- sum(gpt_disgust["covered"], na.rm = TRUE)/nrow(gpt_disgust["covered"])
percent_m_sadness <- sum(gpt_sadness["covered"], na.rm = TRUE)/nrow(gpt_sadness["covered"])
percent_m_enthusiasm <- sum(gpt_enthusiasm["covered"], na.rm = TRUE)/nrow(gpt_enthusiasm["covered"])



library("ggplot2")
library("reshape")
barplots <- data.frame(rows = c("protest", "symbol", "anger", "fear", "disgust", "sadness", "enthusiasm"),  DSL= c(percent_dsl_protest, percent_dsl_symbol, percent_dsl_anger, percent_dsl_fear, percent_dsl_disgust, percent_dsl_sadness, percent_dsl_enthusiasm),   "GPT-4o" = c(percent_m_protest, percent_m_symbol, percent_m_anger, percent_m_fear, percent_m_disgust, percent_m_sadness, percent_m_enthusiasm))
# Barplot
barplots
df2 <- reshape::melt(barplots, id = c("rows"))
df2$variable <- factor(df2$variable, labels = c("DSL", "GPT-4o (SO)"))


ggplot(data = df2, aes(x = rows, y = value)) + facet_wrap(~ variable) + geom_point(alpha = 0.7, size = 3) + labs(title = "Coverage of 95% confidence intervals", y = "Coverage percentage", x = "Independent variables") +theme(plot.title=element_text(face="bold")) + geom_hline(yintercept = 0.95, linetype="dashed", col = "blue") + theme_bw() 






```

```{r}

a <- data.frame("term" = c("protest", "symbol", "anger", "fear", "disgust", "sadness", "enthusiasm"), "estimate" = c(mean(dsl_protest$estimate), mean(dsl_symbol$estimate), mean(dsl_anger$estimate), mean(dsl_fear$estimate), mean(dsl_disgust$estimate), mean(dsl_sadness$estimate), mean(dsl_enthusiasm$estimate)), "std.error" =  c(mean(dsl_protest$std.error), mean(dsl_symbol$std.error), mean(dsl_anger$std.error), mean(dsl_fear$std.error), mean(dsl_disgust$std.error), mean(dsl_sadness$std.error), mean(dsl_enthusiasm$std.error)), "p.value" = c(mean(dsl_protest$p.value), mean(dsl_symbol$p.value), mean(dsl_anger$p.value), mean(dsl_fear$p.value), mean(dsl_disgust$p.value), mean(dsl_sadness$p.value), mean(dsl_enthusiasm$p.value)), "conf.low" = c(mean(dsl_protest$conf.low), mean(dsl_symbol$conf.low), mean(dsl_anger$conf.low), mean(dsl_fear$conf.low), mean(dsl_disgust$conf.low), mean(dsl_sadness$conf.low), mean(dsl_enthusiasm$conf.low)), "conf.high" = c(mean(dsl_protest$conf.high), mean(dsl_symbol$conf.high), mean(dsl_anger$conf.high), mean(dsl_fear$conf.high), mean(dsl_disgust$conf.high), mean(dsl_sadness$conf.high), mean(dsl_enthusiasm$conf.high)))

b <- b[c(10:16),c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high")]

c <- data.frame("term" = c("protest", "symbol", "anger", "fear", "disgust", "sadness", "enthusiasm"), "estimate" = c(mean(gpt_protest$estimate), mean(gpt_symbol$estimate), mean(gpt_anger$estimate), mean(gpt_fear$estimate), mean(gpt_disgust$estimate), mean(gpt_sadness$estimate), mean(gpt_enthusiasm$estimate)), "std.error" =  c(mean(gpt_protest$std.error), mean(gpt_symbol$std.error), mean(gpt_anger$std.error), mean(gpt_fear$std.error), mean(gpt_disgust$std.error), mean(gpt_sadness$std.error), mean(gpt_enthusiasm$std.error)), "p.value" = c(mean(gpt_protest$p.value), mean(gpt_symbol$p.value), mean(gpt_anger$p.value), mean(gpt_fear$p.value), mean(gpt_disgust$p.value), mean(gpt_sadness$p.value), mean(gpt_enthusiasm$p.value)), "conf.low" = c(mean(gpt_protest$conf.low), mean(gpt_symbol$conf.low), mean(gpt_anger$conf.low), mean(gpt_fear$conf.low), mean(gpt_disgust$conf.low), mean(gpt_sadness$conf.low), mean(gpt_enthusiasm$conf.low)), "conf.high" = c(mean(gpt_protest$conf.high), mean(gpt_symbol$conf.high), mean(gpt_anger$conf.high), mean(gpt_fear$conf.high), mean(gpt_disgust$conf.high), mean(gpt_sadness$conf.high), mean(gpt_enthusiasm$conf.high)))



alltogether <- rbind(b, c, a)




alltogether$model <- c("Benchmark", "Benchmark", "Benchmark", "Benchmark", "Benchmark", "Benchmark", "Benchmark","GPT-4o (SO)","GPT-4o (SO)","GPT-4o (SO)", "GPT-4o (SO)", "GPT-4o (SO)", "GPT-4o (SO)", "GPT-4o (SO)", "DSL", "DSL", "DSL", "DSL", "DSL", "DSL", "DSL")

alltogether$model <- factor(alltogether$model, levels = c("DSL", "GPT-4o (SO)", "Benchmark"))

casas1 <- alltogether %>% 
  ggplot(aes(estimate, model, colour = model, shape = model)) + facet_wrap(~ term) +
  scale_shape_manual(values = 0:9) + 
  geom_point(show.legend = FALSE) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), show.legend = FALSE) +
  # add in a dotted line at zero
  geom_vline(xintercept = 0, lty = 2) + coord_cartesian(xlim = c(-0.1, 0.6))+
  labs(title = "Estimate of effect of variable on number of retweets",
       y = NULL) + theme_bw() + scale_color_grey()

```

```{r}

truth_baselines <- data.frame(
  term = c("anger", "disgust", "enthusiasm", "fear", "sadness", "protest", "symbol"),
  truth_value = c(t_anger, t_disgust, t_enthusiasm, t_fear, t_sadness, t_protest, t_symbol)
)




casas2 <- ggplot(data = df2, aes(x = rows, y = value)) + facet_wrap(~ variable) + geom_point(alpha = 0.7, size = 3) + labs(title = "Coverage of 95% confidence intervals", y = "Coverage percentage", x = "Independent variables") +theme(plot.title=element_text(face="bold")) + geom_hline(yintercept = 0.95, linetype="dashed", col = "blue") + theme_bw() +
  theme(
    plot.title = element_text(size = 15), 
    axis.title = element_text(size = 13),               
    axis.text = element_text(size = 13),
    strip.text = element_text(size = 13) # Adjust facet label font size
  ) +
  scale_color_grey()


casas1 <- alltogether %>% 
  ggplot(aes(estimate, model, shape = model)) + facet_wrap(~ term) +
  scale_shape_manual(values = 0:9) + 
  geom_point(show.legend = FALSE) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), show.legend = FALSE) +
  # Add variable-specific dotted lines
  geom_vline(data = truth_baselines, aes(xintercept = truth_value), lty = 2) + coord_cartesian(xlim = c(-0.5, 0.5))+
  labs(title = "Estimate of effect of variable on number of retweets",
       y = NULL, x = "Estimate") + theme_bw() + scale_color_grey() +
  theme(
    plot.title = element_text(size = 15), 
    axis.title = element_text(size = 13),               
    axis.text = element_text(size = 13),
    strip.text = element_text(size = 13) # Adjust facet label font size
  ) +
  scale_color_grey()


require(gridExtra)

grid.arrange(casas1, casas2)
```
